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A model is presented for the characterization of dissipative effects on highly nonlinear waves in one- 
dimensional dry granular media. The model includes three terms: Hertzian, viscoelastic and a term 
proportional to the square of the relative velocity of particles. The model outcomes are confronted 
with different experiments where the granular system is subject to several constraints for different 
materials. Excellent qualitative and quantitative agreement between theory and experiments is 
found. 
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There exist a considerable number of engineering appli- 
cations where surfaces are subjected to contact loading, 
with stress applied over small areas. Nevertheless, it is a 
complex matter to understand the nature of the contact 
of solid surfaces. A first step towards this understand- 
ing started with the work by Hertz [1]. Using potential 
theory, he developed a model for the normal contact on 
non-conforming bodies of elliptical shapes, under several 
assumptions that simplified the problem. One of the hy- 
potheses on which his work was based is that the contact 
between solids is purely elastic. 

A place where Hertz theory finds application is in the 
study of granular matter. This kind of matter is present 
everywhere in Nature and has wide practical importance 
(see e.g., [2]). Hertz theory has been useful in under- 
standing many aspects of granular matter [sl, 
but its applicability is limited because in many situa- 
tions the impact between grains is such that energy dis- 
sipating phenomena become relevant. Such phenomena, 
involving for example elasto-plastic and viscoelastic be- 
havior [5,], are so complex that it is hard to think of a 
closed form force law that may describe them all at once. 

However, energy loss can be measured experimentally, 
at the macroscopic level, by measuring the coefficient of 
restitution, which has been observed to decrease with the 
normal component of the relative impact velocity [6]. In 
this case. Hertz theory fails to reproduce the behavior of 
the coefficient of restitution. In and [H, the authors 
have developed a quasi-static approximation to calculate 
the normal force acting between colliding particles, as- 
suming a viscoelastic force (see also [9] where this force 
law appears in a related context, and [lO| for a first- 
principles derivation of the viscoelastic term). Within 
this approach, where the force acting between beads is a 
combination of Hertz and viscoelastic terms, theory and 
experiment agree for the behavior of the coefficient of 
restitution with velocity [llj. 

In this Letter, a model that combines Hertz theory, a 
viscoelastic force as in Q and [H, and a force propor- 
tional to the square of the relative velocities of beads 
is presented. It is shown that this model reproduces in 
an excellent way the effect of dissipation on a solitary 



wave in stainless steel, brass and polytetraflouroethylene 
(PTFE) obtained by Daraio et al. pjj. It also coincides 
very well with the behavior of solitary wave trains in a 
column of stainless steel beads [l^ and with the descrip- 
tion of incident and reflected solitary waves in a column 
of PTFE balls [isj. To simulate such systems, I follow 
as closely as possible the experimental setup, including 
the reduction by less than 5.5% in mass of beads with 
inserted piezosensors. 

Let Xi{t) represent the displacement of the center of 
the i-th sphere, of mass m^, from its initial equilibrium 
position. The equations of motion that describe the dy- 
namics of N beads, inclined by an angle a, in a gravita- 
tional field are: 

+ 5|6>i_i {sign[Si-i] 5i-if - Oi sign[5i] 5'^^ 
+ mi^sin[a], (1) 

with z = 2,...,7V — 1. As known, the equations of mo- 
tion for the first and last beads differ from Eq. ([T]) 
and are thus not written down here. The notation 
is as follows: Si = Xi — x^+i is the relative velocity 
of beads i and i -\- 1. The overlap between adjacent 
beads is Si = maxjA^^^+i — (x^+i — Xi),0}, ensuring 
that the spheres interact only when in contact. For 
the same reason a step-function in the third term has 
been included, that is, Oi = ^[A^^^+i — (x^+i — Xi)]. 

= {g sin[a] z mi/Ki^i+i)^/^ appears from the pre- 
compression due to the gravitational interaction. The 
expression for the Hertz coupling Kij between beads i 
and j is well known and depends on radii. Young moduli 
and Poisson ratio of beads [5] . Although the form of pa- 
rameter A is known for a binary collision jsj, it is used 
here as a free parameter, like B itself. The set of Eqs. 
([1]) is solved by using an explicit Runge-Kutta method 
of the 5th order with an embedded error estimator, from 
Mathematica. 
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Our explanation for our Ansatz is as follows: after the 
impact, the dynamics becomes a multi-impact problem; 
this produces that the relative velocity of beads i and 
i + 1, may change from < 0, related to an expansion 
phase, to Si > 0, corresponding to a compressional phase. 
The same happens for the relative velocity of beads i — 1 
and i, Si-i. Therefore, by fixing our attention on one 
bead in the chain, say i — th, one has a particle between 
two moving walls (beads i and i -\- 1) and then its dy- 
namics depends on the dynamics of both constraints. It 
is worthwhile to mention that a term proportional to the 
square of the velocity was introduced by Poschl [l5| as 
an attempt to extend the Hertz theory to plastic bodies. 
Also notice that if beads i — 1, i and are all in contact 
at a given instant, one can easily observe that the third 
term can be generically written as 
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2ex'^ -\- K, Xi_^i -\-2r]Xi {xi-i — cr x^+i ) , 



(2) 



where the constants take the values e = 0,1, n = ±1, 
77 = =bl, cr = ±l, depending on the sign of the relative 
velocities 5i and Observe that the force term in- 

cluded in this Letter is a combination of two terms, only 
one of them being dissipative [Ig']. 



In Ref. [12|, the effect of dissipation on the behav- 



ior of solitary waves was clearly shown. Beads made of 
stainless steel, brass and PTFE, and a wall of aluminium 
were used in the experiment. Their Young modulus and 
Poisson ratio are: (i) stainless steel: E = 193 x 10^ Pa; 
V = 0.30; (ii) brass: £; = 115 x 10^ Pa; v = 0.31; (iii) 
PTFE: ^ = 1.46 X 10^ Pa; v = 0.46; and (iv) aluminium: 
^ = 69 X 10^ Pa; u = 0.33. Strikers with the same me- 
chanical properties as beads were used to generate soli- 
tary waves; the force on piezosensors was recorded, and 
the data presented as plots of force as a function of time. 

In Figure 1, the numerical findings from our model are 
compared with those shown in Fig. 1(b) of Ref. [12], that 
is, for a chain composed of = 70 stainless steel beads, 
and impactor velocity vi = 1.77 m/s. Sensors are placed 
in beads 9, 16, 24, 31, 40, 50, 56, and 63. A global time 
shift of 25 jiis of the extracted data was necessary in or- 
der to compare our findings with the experimental data 
from Fig. 1(b) of Ref. [12]. This shift possibly can be as- 
cribed to an experimental time offset. Also, in the data 
received by the author, sensors appear in pairs, say sen- 
sors 2 and 4 (e.g. beads 16 and 31 for stainless steel). 
Thus, in order to compare this experimental data with 
the simulation, the original experimental data has been 
shifted such that the time interval between both max- 
ima in the data is kept fixed, within an error of 0.2 /is. 
Thus, this shift does not represent a change in the time 
of fiight, lying within the experimental errors. Figure 1 
also shows, as an inset, the percentage error between ex- 
periment and simulation. One observes that this error 
is below 20% for the maximum amplitude, getting worse 
when approaching the wings of the signal, which is not 
surprising from the experimental point of view. Figure 
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FIG. 1: The continuous curve shows our numerical findings. 
Open circles and error bars are experimental data extracted 
from Fig. 1(b) of Ref. [12]. The experimental outcome from 
sensors 16, 31, 40 and 56 is also shown. The inset shows 
the percentage error between experiment and simulation for 
beads 16 (triangles) and 40 (circles). 




FIG. 2: The figure shows the force as a function of time for 
beads (a) 16, and (b) 56, for stainless steel. Dashed lines 
correspond to the numerical results. 



2 shows a detail of Figure 1, using original data, for the 
force recorded at beads 16 and 56. As one can see, the 
agreement between simulations and experiment is quite 
impressive. One possible combination of values of the 
parameters, not necessarily the optimal one, that gives a 
quite remarkable coincidence with experiment is A = 800 
and B = -1.9. 
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FIG. 3: Plots showing force as a function of time, for an 
impact velocity of vi = 1.55 m/s, for (a) PTFE (bead 38), 
and (b) brass (bead 49). 



In Figure 3, a striker with initial velocity vi = 1.55 m/s 
impacts on a chain of = 69 PTFE beads and a chain 
of N = 61 brass spheres. Forces are recorded by sensors 
at beads 38 and 49, respectively. The numerical values 
of the parameters that caused theory and experiment to 
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FIG. 4: Scattering of highly nonhnear waves off a wall for 
PTFE beads at positions: (a) 12 and (b) 16. In (c) the sensor 
is at the wall. Dots represent experimental data and the solid 
line the numerical output. 
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FIG. 5: Scattering of highly nonlinear waves off a wall for 
stainless steel beads at positions: (a) 12, and (b) 16. In (c) 
the sensor is at the wall. Data is shown as dots and numerical 
results as a solid line. 



coincide are A = 100 and B = -0.26, for PTFE, and 
A = 770 and B = -2.9, for brass. 

In the following, the simulation findings from our 
model are compared with those from experiments car- 
ried out by Nesterenko et al. with a column of PTFE 
beads [l^ lying on a wall made of brass. There, soli- 
tary waves were generated by impacting the column of 
21 beads, each with radius 2.38 mm and mass 0.123 g, 
with a PTFE ball with the same characteristics and with 
a 2.0 m/s impact velocity; the mechanical properties of 
materials being the same as above. In Figures 4 (a) and 
(b), there appear the incident and reflected perturbations 
recorded at positions 12 and 16, respectively, while Fig- 
ure 4 (c) shows the behavior of the solitary wave at the 
wall. Because of the lack of detailed experimental in- 
formation, global time shifts by 26 /is and 38 /is where 
applied to get coinciding incident perturbations in fig- 
ures (a) and (b), respectively, while for figure (c), the 
shift made was 31 jas. Observe that the simulation fits 
quite well the experimental data. The same numerical 
values for the parameters of the model in case of PTFE 
are used, i.e. A = 100 and B = —0.26. 

To end the tests on the validity of the model, the sim- 
ulation results are compared with the outcomes from ex- 
periments carried out by Nesterenko et al. with a col- 
umn of stainless steel beads [l^. The data, extracted 
from their paper, are used to show that the model re- 
ported here successfully describes the behavior of highly 
nonlinear waves on hard walls. In order to create the non- 
linear waves, a column of 21 stainless steel beads, each 
with a mass 0.45 g, was struck by a cylindric alumina 
impactor, with Young modulus E = 416 x 10^ Pa and 
Poisson ratio u = 0.23. The cylinder has a mass 1.2 g 
and a velocity equal to 0.44 m/s. The force is recorded 
by embedding piezosensor in beads 12 and 16 (measured 



from the top of the column), and at the wall; the piezo- 
gauge at the wall was covered by a brass cover-plate. 
As expected, because the impactor's mass is much larger 
than the mass of beads, trains of highly nonlinear soli- 
tary waves are excited by the impact. Figure 5 shows 
an excellent agreement between model and experiment; 
it must be stressed that the numerical values for the pa- 
rameters are the same as before: A = 800 and B = —1.9. 
In order to compare experiment and simulation, the data 
was shifted 47 /is for figures (a) and (b), and 45 /is for 
figure (c). The small time difference, of the order 7 /is, 
between the pulse amplitudes from the simulation and 
the experiments observed in Figures 4 and 5, (a) and 
(b), can be ascribed to the combined properties of both, 
wall and piezosensor. If a softer wall is assumed, the time 
difference and the amplitude of the reflected wave (and 
the force at the wall) can be reduced. 

Finally, I compare the contribution of the viscoelastic 
and velocity-squared terms, using the force on bead 10. 
The Hertz interaction is by far the leading term and it 
is not included here. In Figure 6, the force on bead 10 
is plotted against time; time units are not made explicit 
because not all of the plots correspond to the same time 
interval. Also, in Figures 6 (a) and 6 (b) force values 
at the r.h.s. correspond to PTFE. In Figure 6 (a), one 
observes that both terms give a contribution of the same 
order, but in some time intervals each acts in opposition 
to the other. A slightly different situation appears in 
case of the reflected perturbations, as seen in Figure 6 
(b). In this case, there is a time interval where both force 
terms reinforce their contribution. Also observe that the 
velocity-squared force term is a continuous, although not 
smooth, function of time. Nevertheless, this is not crucial 
for reproducing the experimental data. 

In conclusion, a model that reproduces experimental 



4 



1.0 

-1.0 

-2.0 
0.15 



-0.15 



//% (a) ^ 










/ \\ 
*_ 


(b) 




, \ 


/ // 



0.1 

0.0 

-0.1 
0.04 



-0.04 



FIG. 6: Figure (a) shows viscoelastic (1) and velocity- squared 
(2) forces as a function of time on bead 10, for beads of stain- 
less steel (dashed (1) and dotted (2) lines), brass (long (1) 
and short (2) dot-dashed) and PTFE (thin (1) and thick (2) 
continuous), with the conditions of Figures 1,2 and 3. Figure 
(b) shows the force terms in case of reflected waves, using the 
data used in Figures 4 and 5, for stainless steel (thin (1) and 
thick (2) continuous) and PTFE (dashed (1) and (2) dotted). 



outcomes for the behavior of highly nonlinear dissipative 
waves, for different materials and under different condi- 
tions, in one-dimensional dry granular media has been 
found. Excellent qualitative and quantitative agreement 
between theory and experiments is found, within the ex- 
perimental errors of measurements of either force ampli- 
tude or time of flight, mostly depending on the Hertz 
interaction. The force term added in this Letter, pro- 
portional to the square of the relative velocity of parti- 
cles, completes a previous physical model composed of 
Hertz and viscoelastic interactions which, without this 
term, is unable to reproduce experimental outcomes for 
the behavior of dissipative highly nonlinear waves. The 
new force term originates from multiple impacts within 
the system and is composed of two parts, one of which 
is dissipative. The model is economical, requiring only 
two parameters that depend on the mechanical proper- 
ties of beads, one of them being well known and with 
clear physical origin. In addition, the parameter values 
do not need to be changed in order to fit simulations with 
experiments carried out under different conditions, and 
are found through a single fit procedure. Of course, if one 
is looking for the optimal set of parameters, one should 
perform a more complete analysis, like the one done in 
[12,]. The model is, of course, not universally valid and, 
for example, is useless for describing the generation of 
oscillatory shock waves in soft materials as PTFE. 
The author is indebted to Prof. C. Daraio for discus- 
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